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Stability of self-similar solutions for gravitational collapse is an important problem to be investi- 
gated from the perspectives of their nature as an attractor, critical phenomena and instability of a 
naked singularity. In this paper we study spherically symmetric non-self-similar perturbations of 
matter and metrics in spherically symmetric self-similar backgrounds. The collapsing matter is as- 
sumed to be a perfect fluid with the equation of state P = ap. We construct a single wave equation 
governing the perturbations, which makes their time evolution in arbitrary self-similar backgrounds 
analytically tractable. Further we propose an analytical application of this master wave equation to 
the stability problem by means of the normal mode analysis for the perturbations having the time 
dependence given by exp (z'wlog \ t\), and present some sufficient conditions for the absence of non- 
' oscillatory unstable normal modes with purely imaginary co. 
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I. INTRODUCTION 

CN 
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Gravitational collapse is a long-standing problem to be investigated in general relativity. Even if the analysis 
• is restricted to a spherically symmetric system, the Einstein equations still remain too complicated to understand 
fully the generic features of the process of gravitational collapse. One may try to overcome such a difficulty by 
considering self-similar spherically symmetric spacetimes, for which the Einstein equations are reduced to a set of 
ordinary differential equations with respect to the one variable z = r/t, where r and f may be comoving radial and 
jy^ , time coordinates. Such self-similar solutions are useful for seeing interesting phenomena in relation to the violation 
of the cosmic censorship due to naked singularity formation and the critical behavior at the threshold of black hole 
\Q formation. Further, as was advertised by Carr (see recent review [1]) as the so-called self-similarity hypothesis, a self- 
similar behavior may be an attractor behavior, that is, it may get dominant near a central dense region as gravitational 
collapse starting from generic initial conditions proceeds to final stages. Though this hypothesis strongly motivates 
us to study extensively self -similar models, its validity should be confirmed through detailed analysis of stability of 
^ ' self-similar behaviors for non-self-similar perturbations. 

One approach to the stability problem may be to study time-evolution of massless test (scalar or electromagnetic) 
fields in self -similar background spacetime [2, 3]. In particular, the possible generation of burst-like emission of 
electromagnetic radiation due to the infinite blue-shift effect at the future Cauchy horizon associated with a central 
naked singularity was discussed in [3] by means of the Green's function technique. Though the test-field approach 
is interesting for examining the instability of the Cauchy horizon associated with a naked singularity, our concern 
in this paper is rather to treat directly non-self-similar perturbations of collapsing matter, which are accompanied 
with metric perturbations. To make our analysis tractable, we consider only spherically symmetric perturbations. 
Further the collapsing matter is assumed to be a perfect fluid with pressure P given by the equation of state P = xp, 
where a is a constant in the range < a < 1 . This is because there exist various interesting classes of self-similar 
spherically symmetric perfect fluid solutions (see [4] and [5] for example), and their stability has great implications 
to astrophysical problems. 

In perfect fluid dynamics a sonic point located at z = z s < 0, whose value is dependent on the equation of state 
parameter a, plays the role of a characteristic surface, instead of the null Cauchy surface associated with the central 
singularity for massless test fields. Physically allowable self-similar solutions are constrained by regularity (namely, 
at least C 1 ) at the sonic point, and the family of transonic solutions given in the subsonic region between the regular 
center z = 0~ and the sonic point z = z s contains one parameter D in addition to the parameter a. The classification 
of the transonic solutions as well as the clear presentation of the global spacetime structure has been done [6, 7]. 

The banded structure of the parameter D is a notable point of such a classification by the total number N of 
recurrent bounces and recollapses which may occur subsequently to an initial collapse. The allowed range of D for 
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the N-th class of solutions is limited to a finite band D2N < D < D2N+1 (N = 0, 1, 2, • • • ). If the solutions are required 
to be analytic (namely, C°°) at the sonic point, the parameter D is further restricted to some discrete values lying in 
the N-th band. For example, there exist only two no-bounce solutions (corresponding to N = 0) analytic at the sonic 
point [6]. They are called the flat Friedmann solution and the general relativistic Larson-Penston solution, which 
describe a homogeneous collapse and an inhomogeneous collapse, respectively It is also remarkable that the perfect 
fluid critical solution [8, 9] corresponding to the threshold of black hole formation is the N = 1 solution analytic at 
the sonic point. 

Several works have been devoted to numerical analysis of the stability problem for the self-similar perfect fluid 
solutions analytic at the sonic point. In particular, it was numerically confirmed for the parameter range < a < 
0.036 that the flat Friedmann solution and the general relativistic Larson-Penston solution can act as an attractor in 
generic gravitational collapse [10], while the critical solution for the parameter range < a < 0.89 was numerically 
shown to be unstable [11, 12, 13]. (It is interesting to note that the general relativistic Larson-Penston solution for the 
parameter range < a < 0.0105 describes naked singularity formation [6], and thus it is a serious counterexample 
against the cosmic censorship.) Some of the N > 2 solutions were also numerically found to be unstable f or < a < 
0.036 [10] and a = 1/3 [13]. These numerical results obtained for limited values of a and D suggest that only the 
N = solutions can remain stable. Such an interrelation between the stability and the occurrence of bounces is also 
an interesting issue to be further examined in the light of the analytical treatment of the stability problem for any a 
and any allowed range of D. 

We start from a brief review of general spherically symmetric perfect fluid system in Sec. II and mention some key 
properties of the self-similar solutions in Sec. III. Further, in Sec. IV, the set of several linear equations for non-self- 
similar matter and metric perturbations is derived, and the perturbation equations are reduced to a single second- 
order wave equation. This construction of the single wave equation governing the matter and metric perturbations 
will allow us to study analytically time evolution of perturbations in a well-developed standard manner, for instance, 
through the Green's function technique. As a preliminary application of this master wave equation to the stability 
problem, in Sec. V, we try to solve the eigenvalue problem of the spectral parameter co for perturbations having the 
time dependence exp (ico log |f |) and satisfying the boundary conditions (which exclude the kink modes studied in 
[14]) at the regular center z = 0~ and the sonic point z = z s . Then, as the main result in Sec. V, some sufficient 
conditions for the absence of unstable modes with purely imaginary w in self -similar backgrounds (non-bouncing at 
least in the subsonic region) are presented. In the final section, we summarize the results obtained in this paper and 
discuss the consistency of our analytical scheme with the numerical results obtained in previous studies. 



II. SPHERICALLY SYMMETRIC PERFECT FLUID SYSTEM 



In this section we briefly review the general spherically symmetric perfect fluid system as a first step to study 
matter and metric perturbations in self -similar backgrounds, mainly referring to the paper [10]. The spherically 
symmetric line element is given by 

ds 2 = -e 2v ^dt 2 + e 2A ^dr 2 + R 2 (t,r) (dQ 2 + sin 2 Odcp 2 ) , (2.1) 

with the comoving coordinates t and r. The collapsing matter is a perfect fluid whose energy momentum tensor is 
given by 

T ab = (p + P)u a u b + Pg ah , (2.2) 

where p, P and the vector u a are energy density, pressure and fluid four velocity, respectively. As was mentioned in 
Sec. I, we assume that its equation of state is 

P = ap , (2.3) 

using a constant a lying in the range < a. < 1. To discuss self -similar behaviors in the following, we use a new 
variable z defined by z = r/t, instead of r. In addition we also introduce the following dimensionless functions: 

rj(t,r) = 8nr 2 p, (2.4) 
S(t,r) = j. (2.5) 
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From the Einstein's field equations, we can obtain the four equations governing the functions v, A, S and rj. By 
virtue of the choice of the comoving coordinates, the two equations lead to the relations 

e v = C v (t)(z 2 ) a/(1+a V a/(1+a) , (2.6) 
e A = C x {r)rj- l/{ - l+0 ^S- 2 , (2.7) 

where C v and C A are arbitrary functions. Thus the remaining two equations become the equations for only the two 
unknown functions S and rj and are written by 

M + M' = rjS 2 {S + S') , (2.8) 
M- M' = -<xriS 2 (S- S') , (2.9) 

where the dot and the prime represent the partial derivative with respect to log |t| and log |z|, respectively. The 
function M introduced in Eqs. (2.8) and (2.9) is defined by 

M(f,z) = S {l + e~ 2v z 2 (S - S') 2 - e~ 2X (S + S') 2 } . (2.10) 

Note that the function M is the dimensionless function related to the Misner-Sharp mass m as M — 2m/r. In the 
following Eqs. (2.8) and (2.9) will act as the basic equations to derive the perturbation equations. 
Finally in this section we define another important quantity V as 

V(t,z) = ze A ~ v . (2.11) 
This is interpreted as the velocity of a z = const surface relative to the fluid element. 



III. BACKGROUND SELF-SIMILAR SPACETIMES 



Now we point out some features of spherically symmetric self-similar backgrounds. In particular we focus our 
concern on their asymptotic behaviors near the regular center and near the sonic point, at which some boundary 
conditions are imposed on the perturbations. The self-similarity considered here means that 

v = v{z), A = A(z), S = S(z), rj = rj(z), (3.1) 

which require C v and C A be constant. Then, Eqs. (2.8), (2.9) and (2.10) are reduced to the ordinary differential 
equations, which are written as 

M + M' = r/S 2 (S + S') , (3.2) 
M' = -arjS 2 S' , (3.3) 
M(z) = S {l + e- 2v z 2 S' 2 - e~ 2A (S + S') 2 } . (3.4) 

It can be easily found that the Ricci scalar constructed by the self-similar metrics indefinitely increases as t ap- 
proaches zero along a z = const line. Thus, in the limit t — > 0, a singularity will appear at the center r = in 
the self-similar spacetime. Hereafter we are interested in the time evolution of the system before the singularity 
formation, i.e., t < (i.e., z < 0). 

Because the derivative of the function S along a r = const line leads to 

d S dS 

z Tz ^ ^It' (3 - 5) 

the inequality S' > (or S' < 0) means a local expansion (or local contraction) of a fluid shell. Taking account of such 
a physical meaning of the function S' , we here introduce our original expression for Eqs. (3.2), (3.3) and (3.4) which 
is useful for seeing the motion of fluid; we rewrite Eqs. (3.2), (3.3) and (3.4) in terms of the function S' (or A' = S' /S) 
and the velocity V, instead of the functions rj and S. From Eqs. (2.6) and (2.7) we find that the function V is related 
with the functions rj and S by 

V = -C A C- 1 (-z)( 1 -") / ( 1+a )^' l - 1 ) / ( 1 + a )s- 2 . (3.6) 
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In addition the function M can be written in terms of the functions A' and V as 

M = S [l + x{v 2 A /2 - (l+A') 2 }] , (3.7) 

where the function X is defined by 

X = S 2 e~ 2A . (3.8) 
Because of Eq. (2.7), the function X is also rewritten as 

X - C-y /(^'S 6 . (3.9) 
Using Eqs. (3.6), (3.7) and (3.9), we can reduce Eqs. (3.2), (3.3) and (3.4) to 

V (1 - a 2 ) A" + (1 + a) (1 - 5a) A' 2 - 2(a 2 + 2a - I) A' + 1 - a 



V (1 + «){1 + (1 + «)A'} 

X' _ -2(l + a)A" + 6a(l + a)A /2 +4aA' 



(3.10) 
(3.11) 



X (l + a){l + (l + a)A'} 

2X(1 + a) (y 2 - aj A" + (1 - a)XV 2 A' {3(1 + a) A' + 2} 

-X(l + A') {(1 + «)(1 + 7oc)A' + a 2 + 6a + l} + (l + a) 2 = 0. (3.12) 

Note that the use of the functions X and A' instead of M and S' prevents the functions S and ^ and the variable 
z from explicitly appearing in Eqs. (3.10), (3.11) and (3.12) to reduce the field equations to a simpler autonomous 
system. In this paper we will mainly represent the self-similar solutions by the functions V, A' and X. 

One of the requirements for the solutions to describe gravitational collapse beginning with the regular initial data 
is their regularity at the center r = during t < (i.e., z = 0~). The approximate behaviors of the solutions S, rj 
and M satisfying the regularity condition near the center z = 0" are given in [10], from which the behavior of the 
function V near the regular center is found to be 

V ~ -(1 - p) 2/3 d[ /3 (2D)-P /2 (-z) 1 -P , (3.13) 
where p and D are constants defined as 

' S 3(lW (314) 

D = 4np(t,0)t 2 . (3.15) 

In addition we can also find that at the regular center 

A' = -p, X=(l-p)~ 2 . (3.16) 

Because the constant C\ is just the freedom of rescaling the radial coordinate, Eq. (3.13) means that the self-similar 
solutions with the regular center are characterized by one parameter D for a given a. 

From Eq. (3.12) we see that the self-similar perfect fluid system may become singular when the velocity of a 
z = const surface relative to the fluid is equal to the sound speed, i.e., V 2 = a.. We call this singular point in the 
equation a sonic point and denote its location by z = z s . For convenience we hereafter denote the value of a function 
at the sonic point by attaching the suffix "S" to the function (e.g., A'(z s ) = A' s ). 

Because of Eq. (3.13), the function V becomes zero at z = 0~. This means that the flow speed is subsonic at 
t = -oo. However, with the lapse of time, the flow speed of a collapsing shell with a constant r will get larger. Thus 
in gravitational collapse the velocity V (or A' and X) can smoothly decrease to a supersonic value V < — *fk~ beyond 
the sonic point z = z s as z decreases. The asymptotic behaviors of the self-similar solutions near the sonic point and 
the transonic conditions were examined in e.g., [15] and [6]. We here briefly review such studies, using the result 
obtained in Appendix A, in which the asymptotic behaviors of our original set of solutions V, A 1 and X near the 
sonic point are examined. 

The transonic solutions, namely, the solutions which are at least smooth at the sonic point, are well parametrized 
by A' s (or V' s ) in the region near the sonic point. Note that V' s is necessarily negative. In addition the terms propor- 
tional to (£ — £ s )* with non-integer power index % generally appear in the expansion of the solutions near the sonic 
point, where £ is defined by £ = log(-z), and the power index % is given by 

X = -1 - # (3-17) 
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In order for the solutions to be smooth at the sonic point, the power index x must be larger than unity, otherwise the 
coefficients (V x , AL and X x written in Appendix A) of the terms including (£ — £ s ) x must vanish by a fine tuning. 
For the former condition the solutions are smooth but not analytic at the sonic point, while the solutions are analytic 
at the sonic point for the latter condition. 

From Eq. (3.17) we see that the transonic condition % > 1 restrict the allowed range of Vg. Because for a given a. 
the value of V' s and the coefficients of the terms including (£ — £ s )* must depend on the parameter D, the transonic 
conditions seem to also restrict the range of possible values of D. In fact, as was mentioned in Sec. I, the banded 
structure D 2 n < D < D 2 n+i {N = 0, 1,2, • • • ) of the allowed range of D for the transonic solutions and the dis- 
cretization of the allowed values of D for the solutions analytic at the sonic point were numerically found. One of 
such discrete values of D is Dp = 2/3(1 + a) 2 , for which we can find the solution for Eqs. (3.10), (3.11) and (3.12) 
written by 

V = V F (-z) 1 ~P / A' = -p, X=(l-p) 2 (3.18) 

with the constant Vp dependent only on a, which is called the flat Friedmann solution. It should be noted that Dp 
lies in the 0-th band, namely, in the range Do < Dp < D\. 

As was also mentioned in Sec. I, it is remarkable that the number of zeros of the solution A' , namely, the total 
number of turns of radial motion of a fluid shell from a contraction to an expansion and from an expansion to a 
contraction, for D lying in the N-th band is equal to N. 1 For example, as seen from Eq. (3.18), the function A' for the 
flat Friedmann solution does not become zero at any z, namely, the flat Friedmann solution describes collapse with 
no-bounce. 



IV. MASTER WAVE EQUATION GOVERNING PERTURBATIONS 

Now we consider spherically symmetric non-self -similar perturbations in self-similar backgrounds describing the 
gravitational collapse of a perfect fluid. Let us begin with expressing the solutions of the field equations (2.8) and 
(2.9) as 

S(f,z) = S B (z){l + eSi(r,z) + 0(e 2 )} , rj(t,z) = rj B (z) {l + e m (t,z) + 0(e 2 )} , 

M(t,z) = M B (z){l + eMi(r,z) + 0(e 2 )} (4.1) 

with a small parameter e. We use the freedom of the coordinate transformation for t and r to require the constants 
C v and C\ in Eqs. (2.6) and (2.7) be non-perturbed, namely, 

C v = C vB , C A = C AB . (4.2) 

Then the perturbations for the metrics v and A are written by Si and r\\ via Eqs. (2.6) and (2.7). In addition we can 
easily write the perturbation r\\ by Si and Mi, using the perturbation equations (up to the linear order of e) derived 
from Eqs. (2.8), (2.9) and (2.10) which do not contain the derivative of rj. We thus obtain the following two first-order 
partial differential equations for Si and M\: 

Pi(z)Mi + P 2 (z)Mi + P 3 (z)Si + P 4 (z)Si + P 5 (z)Si = 0, (4.3) 
Pi(z)Mi - P 1 (z)M[ + Q 2 (z)Mi + Q 3 (z)Si + Q 4 (z)Si + Q 5 (z)S 1 = (4.4) 

with the coefficients P; and Q; written by the background self-similar solutions, which are given in Appendix B. In 
this paper we focus our analysis on the perturbations Si and Mi. Hereafter we omit the suffix "B" attached to the 
background self -similar solutions for simplicity. 

The complicated form of Eqs. (4.3) and (4.4) may enforce us to study numerically time evolution of the perturba- 
tions. Nevertheless, we can expect that it becomes easier to understand analytically its essential features if the two 



1 In [6] Ori and Piran originally revealed the relation between the number of the zeros of the radial velocity u r of the fluid in the non-comoving 
coordinate system and the rank of the permitted band which the solution belongs to. The velocity u r can be related to the quantities in the 
comoving coordinate system used in this paper as u r = VX^-^A' (see [6] for the coordinate transformation from the non-comoving coordinate 
system to the comoving one). Because the functions V and X cannot be zero for z < by definition, the zeros of the velocity u r correspond to 
those of the function A'. 
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equations are reduced to a single second-order wave equation. Here we adopt such a strategy for our perturbation 
analysis. For this purpose we introduce the function Y written as 

Y(t,z) = S 1 (t,z)-f(z)M 1 (t,z) (4.5) 

with an arbitrary function f(z). Then, from Eqs. (4.3) and (4.4) we can obtain the two first-order partial differential 
equations for Y and Mi with the coefficients containing the function /. It is easy to see that if the function / is given 
by 

f _ f {l + (l + «)Ai}{^VA' T (l + A>)} 

1 f± ~ y/£{VA'±y/£(l + A')} ' K ' 

from the equations for Y and Mi we obtain the relation in which both M' and M are eliminated as follows, 

Y(z)M x = Ti(z)Y + T 2 (z)Y' + T 3 (z)Y, (4.7) 

where the coefficients Y and T, are written by the functions P,-, Q, and / (see Appendix B). By using Eq. (4.7) to 
eliminate Mi from the equations for Y and Mi, we arrive at the following single equation for Y only: 

Y-2Y'+ (l-^) Y" + _Ri(z)Y + R 2 (z)Y' + R 3 (z)Y = 0, (4.8) 

where the coefficients are also written by the functions P,- and Q, and /. As will be shown in Sec. V B, a gauge mode 
corresponding to an infinitesimal change of the background self-similar solution due to a transformation of the time 
coordinate t is allowed as a solution of Eq. (4.8). Because no other choice of / different from Eq. (4.6) is possible for 
obtaining the second-order wave equation (4.8), our scheme is not equivalent to a so-called gauge-invariant approach 
for perturbations. 

Here we consider the transformations of the variables in Eq. (4.8), which are given by t — ► u = log(— t) + /(£) and 
where the functions / and £ are defined as 

The new spacial variable £ runs from zero to positive infinity in the region between the regular center £ — > — oo and 
the sonic point £ —> £ s . By virtue of the transformation (4.9), Eq. (4.8) is reduced to the standard form of the wave 
equation: 

Y MH - Y K + W(£)Y, „ + F(£)Y c + !I(OY = . (4.10) 

We can explicitly write the functions W, F and U by the background self-similar solutions V, A 1 and X, which are 
shown in Appendix B. 

Though in general we cannot see the detailed dependence of the functions W, F and U on the variable £ without 
numerically solving the field equations (3.10), (3.11) and (3.12), we can analytically find their general behaviors near 
the boundaries. In the following calculations we choose the function / as /+ because by virtue of this choice the 
behaviors of the functions W, F and U near the sonic point (at which V = - s/oc) become much simpler. Because of 
Eq. (3.13) and the definition of the variable £, the asymptotic behavior of the function V near the regular center in 
terms of the variable Z, is 

V.f^f. (4.H) 
3(1 + ol) 

Using this relation, we find that the functions W, F and LI become infinitely large near the regular center as 

W ~ X -, (4.12) 

F - ~\> ( 4 - 13 ) 
U * |r, (4.14) 
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FIG. 1: Variations of W, F and U for the flat Friedmann solution corresponding to ct = 0.01 (solid line), a = 1/3 (broken line) and 
a = 1 (dotted line). 



without depending on the parameters D and a. In addition, from Eqs. (A7), (A8) and (A10), we can find that the 
values W s , F s and U s of the functions W, F and LI at the sonic point are given by 

W s = -F s = 2 ( 1 - ^L) , (4.15) 



Us = 0. (4.16) 

Taking account of such behaviors of the functions W, F and U in Eq. (4.10) near the boundaries, we will consider the 
boundary conditions for Y in the next section. 

In addition we should mention the behaviors of the functions W, F and U near a point at which A' = 0, namely, 
radial motion of matter stops. Let us denote this point by £ = £ - The functions W, F and U apparently diverge at 
£ = £o m proportion to (£ - £o) _1 - m particular, by using the field equation (3.12) to see the leading term of A' near 
£ = £o/ we nn£ i mat the leading form of the function F near £ = £ is 

(4.17) 



£-Co ' 

Then, it is easy to check that the solution Y of Eq. (4.10) should remain finite or vanish in proportion to (£ - £o) 2 m 
the limit £ — ► £ - Therefore the wave equation (4.10) is applicable even to a background self-similar solution for the 
parameter D lying in the permitted band labelled by N > 1, namely, for a class in which at least one bounce of a 
fluid shell in gravitational collapse is allowed. 

In the case that the background self-similar solution is the flat Friedmann solution, which is given by Eq. (3.18), 
we can explicitly write the functions W, F and U as functions of z. In Fig. 1 we draw their variations as functions of 
x defined by 

* = (4-18) 

which becomes zero at the regular center and unity at the sonic point. The functions W and — F are monotonic and 
always positive for any a. in the range < a < 1. In addition the function U is monotonic and always positive for 
1/9 < a < 1, but not for < a. < 1/9, which will be understood by Eqs. (4.14), (4.16) and the derivative 

2(l + 3«)(l-9«) 

UA1) - — — ' ( ] 
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though the range of x for which U < is too small to be seen in Fig. 1. We expect that such variations of the 
functions W, F and U shown in Fig. 1 may be typical for self-similar solutions for the parameter D lying in the 0-th 
band, namely for solutions describing collapse with no-bounce. 

It should be emphasized that the wave equation (4.10) is applicable to spherical perturbations in any spherically 
symmetric self-similar perfect fluid solutions (in particular, solutions which is smooth but may not be analytic at the 
sonic point) describing the gravitational collapse. Therefore we expect the wave equation to be useful for revealing 
time evolution of the perturbations for various types of models of self-similar collapse covering a wide range of the 
parameter values D and a. As an important application of the wave equation (4.10), we will consider the stability 
problem for self -similar collapse in the next section. 

V. NORMAL MODE ANALYSIS 

As was mentioned in Sec. I, the stability of self-similar backgrounds for the perturbations has been numerically 
studied in many works. Here we would like to use the wave equation (4.10) for an analytical approach to the stability 
problem. Though our analysis done in this section remains preliminary for future investigations, some clear proofs 
supporting the numerical results will be successfully presented. 

A. boundary conditions 

By virtue of the self-similarity of the background, (namely, the fact that the functions W, F and U depend on z 
only) we can reduce the wave equation (4.10) to the ordinary differential equation as follows, 

tptf - Fip£ + (w 2 - icoW -U S jip = 0, (5.1) 

by assuming that the wave function Y is written as 

Y(u,Q = ip(C,a>)j*> u . (5.2) 

Owing to the behaviors of W, F and U mentioned in the previous section, the general solution (or its higher deriva- 
tive) of Eq. (5.1) may become singular at the regular center or at the sonic point. Therefore we study the conditions 
which the function xf> should satisfy at the boundaries. 

From Eqs. (4.12), (4.13) and (4.14), we obtain the asymptotic behavior of the general solution of Eq. (5.1) near the 
regular center as 

^QM^ + QMr 3 • (5-3) 
The divergence of ip implicates that of M\ because of Eq. (4.7). Thus we require that ip becomes 

(5.4) 

near the center. Note that this boundary condition at the regular center is valid for any self-similar backgrounds. 
Because of Eqs. (4.15) and (4.16), near the sonic point the general solution of Eq. (5.1) behaves as 

ip ~ C 3 (cv)e iw Z + C 4 (co)e-^ +i ^ . (5.5) 

Because near the sonic point we have 

(5.6) 

Eq. (5.5) leads to 

Y ~ {c 3 {u>) + C i (co)e-^ +2i ^ e ^ lo §(- f ) , (5.7) 

from which we note that unstable modes (i.e., modes which grow as t — > 0) correspond to solutions Y for co whose 
imaginary part is positive. The boundary condition imposed on Y at the sonic point may be a subtle problem. To 
see this, let us consider the leading behavior of perturbed physical quantities (for instance, near the sonic point. 
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Using Eq. (4.7) and the perturbation equation obtained from Eq. (2.8), we find that r}\ contains the second derivative 
Y". Then, by taking account of Eqs. (5.7), (3.17) and (4.15) we can express the leading behavior of y}\ as 

Vl (t,z) * {riu(a>)+Vib(o>) - ^~^ /v ^} e^(-t) , (5 . 8) 

which should be compared with the background solution r/g- If 7b is analytic at the sonic point, it is plausible that 
the perturbation r\\ is also required to be analytic there, that is, r\\\, is required to be zero. However, the background 
solution may be regular (i.e., C 1 ), but not be analytic there. Then, the non-analytic term in proportion to (£ — £ s )* 
appears in rj B . This fact and Eq. (5.8) mean that for unstable modes (i.e., modes with Im(o;) > 0), the perturbation 
becomes less regular than the background quantity because V s ' is negative. One may claim that the higher derivatives 
of 7/i must remain smaller than the corresponding higher derivatives of f/g (like r\\ for f/g analytic at the sonic point). 
Then, rj lb in Eq. (5.8) should be zero even if the background is not analytic at the sonic point, and the solution for 
Eq. (5.1) must be derived under the condition that the leading behavior of \p near the sonic point is given by 

xp ~ e iw ^ . (5.9) 

Though whether the boundary condition (5.9) is physically necessary is not so conclusive, hereafter we use it as well 
as the boundary condition Eq. (5.4) to set up the eigenvalue problem for Eq. (5.1) and to obtain the normal modes 
and the discrete eigenvalues denoted by ip n and co„. 



B. gauge mode 

We can find an exact solution xpg for Eq. (5.1) written by 



if the spectral parameter to is equal to 



_ c(l+a)VA* j 

VA' + ^(1 + A'f ' (5 ' 1U) 



CO = LO s = - 1 , (5.11) 



where c is an arbitrary constant. We can easily find that the solution xpg satisfies the boundary conditions. Hence co g 
is one of the eigenvalues. This normal mode corresponds to the correction of Sg and Mg generated by the following 
infinitesimal coordinate transformation: 

(-t) - (-f) {l + ec{-tf**} , r -+ r , (5.12) 

which means 

(_ 2 )^(-z){l_ ec (_f )''<"*} . (5.13) 

Therefore we call the solution ip g a gauge mode. Though the gauge mode is obviously unphysical, the presence of the 
solution explicitly written by the background solutions V and A' will become mathematically useful for analysing 
Eq. (5.1), as will be seen later. 



C. Constraints for unstable modes 



In this subsection we present some analyses of Eq. (5.1), which is slightly different from the well-studied Strum- 
Liouville type differential equation. Because we are particularly interested in the stability of the self -similar solutions, 
we examine the possibility of the existence of unstable normal modes. Hereafter we focus our investigation on the 
eigenvalues for the background solutions in which the function A' does not become zero at least in the range between 
the regular center and the sonic point, to avoid the mathematical complexity caused by the singular behavior of the 
functions W, F and U at A' = 0. Such solutions include all the solutions for the parameter D lying in the 0-th band, 
namely, for the class in which no bouncing motion of a fluid shell occurs in collapse. In addition the critical solution 
for 0.61 < a < 1 is also included, because the zero of the function A' is only in the supersonic region [16]. Therefore, 
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the analysis done in this subsection may be applicable to various types of solutions describing collapse accompanied 
with some bounces and recollapses. 

Now let us begin with the introduction of the function tp defined by 

4>(Z,a>) = il>e- iu Z (5.14) 

instead of tp, because normal modes denoted by <p n = xp„ exp {—icu n Q become constant at the sonic point. In addition 
it will be mathematically convenient to use the variable x instead of £ in the following calculations, which can cover 
the whole region between the regular center and the sonic point in the finite range < x < 1, irrespective of the 
parameter values a and D. Using the function <p n and the variable x, we rewrite Eq. (5.1) as 

y/*x L. r V'{x 1 + 1) V"(l-x 2 )l 
<pn,xx - 2T \ llC0n - p+ V 2 + — 17/ " f ^ 

V'(l-X z ) I ^JOLX 2 - V'x J 

2 

OCX 

- y/2 (1 _ x 2)2 {i"n(F + W) + U}<p n =0. (5.15) 

Moreover, denoting the real part and the imaginary part of the eigenvalues to n , respectively by f>„ and y n , we 
rewrite Eq. (5.15) to the following equation: 

( h <Pn,x), x yfif^ Hn* ^(^2)2 ^ (F + W) + U} h<p n = , (5.16) 
where the function h is given by 

** ^" f. 27 „_ f + zv^ + ni^)\ (5 . 17) 



fc - V'il-x 2 ) \ ^LX 2 V'x 

and can be always positive in the range < x < 1. Following the standard procedure, we consider to multiply 
Eq. (5.16) by the complex conjugate <p* of <p n and obtain 

<P*n ( h 4>n,x),x + <Pn {H*n,x) iX ~ y,(i^ x 2) h ( ( Pn<Pn,x ~ <pn<p* n ,x) 

O rv 

{U- 7n (F + W)}h\cp n \ 2 = 0. (5.18) 



V' 2 (l-x 2 ) 2 

Let us integrate Eq. (5.18), taking account of the boundary condition such that <p„ <~ x near the regular center x = 0. 
Because the asymptotic behavior of the function h near the regular center is 



~ x 3 , (5.19) 



we have 



dx 



h(\<Pnf) x = h \\<p n , x \ 2 + ^* x2)2 {U- 7n (F + W)}\cp n \ 

+2ifin Jo W0xtj k ^ n<pn ' x dx ■ (5 ' 20) 

We now apply the boundary condition at the sonic point to Eq. (5.20). Using Eq. (4.15), we find that the asymptotic 
behavior of the function h near the sonic point is 

h~(l-xy a , (5.21) 

where a is a constant defined as 

«=^( 7n -l). (5.22) 
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For 7n > 1 (i.e., a < because V s ' is negative), the integrations in Eq. (5.20) remain finite at x = 1, and the left hand 
side of Eq. (5.20) vanishes there, because <p n and its derivative become constant there. Thus we obtain the relation 



J 



2 

r, (XT O 

I'M" + y,2 {l _ x 2 ) 2 i U -T»( F + W )i W 



dx 



+//5 " jf (^"<* ~ fati*) dx = ( 5 - 23 ) 

for 7„ > 1. This clearly shows that if the inequality 

U-<y n (F + W) >0 (5.24) 

holds for 7„ > 1 and any x in the range < x < 1, then there exists no normal mode with y n > 1 and /3„ = 0. 
It is noteworthy that for j n > 1, Eq. (5.24) is satisfied for the flat Friedmann solution given by any a. in the range 
< oc < 1. 

Though we cannot derive a condition for the absence of unstable normal modes with j n > 1 and /3„ 7^ from 
Eq. (5.23), we can derive another constraint for normal modes <p„ with j6 n 7^ from the multiplication of <p* to 
Eq. (5.16) as follows, 

<Pn {Hn, X ), x - <Pn (Wn,,) , x - " y^W ^ + = " ' ^ 

For normal modes with j n > 1, the integration of Eq. (5.25) gives 

Using Eq. (5.26), we can estimate the second term of the left hand side of Eq. (5.23). Then, the condition (5.23) for 
normal modes with j n > 1 may be shown to be incompatible with the additional condition (5.26), for instance, if the 
background is the flat Friedmann solution. This is an interesting problem to be fully examined in future works. 

Further, we note that if oscillatory unstable normal modes (i.e., unstable normal modes with j6„ 7^ 0) with 7 n > 
W s /2 (i.e., a < —1) exist, Eq. (5.26) allows an estimation of y n . In fact, the partial integration of the term including 
(\<pn\ 2 ),x in Eq. (5.26) becomes possible for j n > W s /2, and we obtain 



[W^xj^ 2dx _ 



r z . . .o . 2 



This means that -y n should be equal to a half of an average W av of the function W. 

We have presented Eq. (5.24) as a sufficient condition for the absence of non-oscillatory unstable modes (i.e., 
unstable modes with f>„ = 0) with j n > 1. Special attention to the case /5 n = will be meaningful because all the 
unstable normal modes numerically found in previous works were non-oscillatory [10, 11, 13]. Therefore the next 
task should be to analyze normal modes with j6„ = and < 7« < 1. It should be noted that there exists a gauge 
mode <p g = ipgexp (-iLOgQ as an unphysical normal mode with f> n = and < 7„ < 1. Therefore, to obtain the 
conditions for the physical normal modes tp n , we must consider a difference between <p n and <p s . For this purpose, 
we multiply Eq. (5.15) for <p n by <p s and Eq. (5.15) for <p s by <p n and subtract the obtained two equations. Then we 
have 

- + ( y_?>jf + ^B^} = " ■ <*» 

where -y g is the imaginary part of co g , which is given by Eq. (5.11). By virtue of the boundary condition at the regular 
center, the integration of Eq. (5.28) leads to 

* (Mn, - = '(in ~ 7 S ) J* { ^ + % F _ + X 7 ) ) ] ****** ■ (5-29) 
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Now let us turn our attention to the asymptotic behaviors of the both sides of Eq. (5.29) near x = 1. Because of 
Eq. (5.21), the left hand side of Eq. (5.29) near x = 1 is given by 

h (<p g <p n , x - <p n <p g/X ) = C 5 M(1 - x)~ a +0 ((1 - x)" fl+1 ) . (5.30) 

If 7„ is in the range 1 + {Vg/y/oi) < J n < 1 (i-e., < fl < 1), the subleading term proportional to (1 - x)~ a+1 must 
vanish in the limit x — > 1. This is a key point to discuss the absence of non-oscillatory normal modes with 7„ giving 
the range < a < 1 by comparing Eq. (5.30) with the right hand side of Eq. (5.29). 
Let us introduce the function H defined as 

= ^(1-xf f 2fex y^x(F + W) 1 

y(i + x) \ ^ + v>(i-x*) j h<pg ' {5Jl) 

which is the function defined by the background solution. It becomes zero at x = and a finite value (denoted by 
H s ) at x = 1. Using the function H, we can rewrite the integral in the right hand side of Eq. (5.29) into the form 

[* VEx f 2<p glX ^x(F + W) \ 

= ^(1 - *)">„ - ^ jf*(l - x)-"- 1 1(1 - x)tf>„,* - « - 1^„| ix . (5.32) 

For < a < 1, the integral in the second term of the right hand side of Eq. (5.32) must vanish in the limit x — > 1 for 
consistency between Eq. (5.30) and Eq. (5.32). Thus we obtain the condition for non-oscillatory normal modes with 
7n giving < a < 1 as follows, 



Jo 



H 



^ (1 - x)- a ~ L I (1 - x)<p n , x -^--lj^* = 0. (5.33) 

It is clear that the condition (5.33) cannot be satisfied if H/H s < 1 and <p n ,x/<pn > for any x in the range < x < 1. 
Fortunately, we can show that if Eq. (5.24) is satisfied, the positivity of <p n ,x/(pn is assured. To see this, we firstly note 
that Eq. (5.15) requires that if (p n becomes locally maximum at m points x = x t (satisfying < X\ < ■ ■ ■ < x m < 1), the 
ratio (p n ,xx/<pn must be positive there. However, recalling that (p n ~ x near x = 0, we can claim that the ratio <p n ,xx/ (pn 
must be negative at x — X\. This contradiction is caused by the assumption of the existence of the points at which 
<pn,x = 0. Hence we can conclude that if a background self-similar system satisfies both the conditions H/H s < 1 
and (5.24) for any x in the range < x < 1, there exists no normal mode with j n giving < a < 1. 

It is interesting to check whether the flat Friedmann solution can satisfy such conditions. Because the ratio H /X /H 
is given by 



H, x _ 6(1 - a)x 2 + 10(1 + 3a)x + 3(1 + 3a) 2 - 3(1 + cc)x{l + 3oc + 2x)<y n 
IT ~ (l+3a)x(l + x)(l + 3« + 2x) 



(5.34) 



it is easy to understand that if 7„ < 1, the function H has no extremum in the range < x < 1, and hence we 
obtain H/H s < 1 there. Further we can find that for 7„ > the second derivative of U — j n (F + W) with respect to 
x become always positive in the subsonic region. Hence, the positivity of U — j„ (F + W) in the subsonic region is 
assured if its derivative with respect to x is negative at x = 1, which leads to 

(l+3«)(l-9«) „.„ nx 

7n> q= — ~A: t 1 ■ (5.35) 

' v 9a 2 + 18a + l v ; 

Note that the value of 1 + (Vg/^/a) for the flat Friedmann solution is equal to p, which is given in Eq. (3.14), and 
that the equality p = q holds at a = a c w 0.01879. Thus the absence of non-oscillatory normal modes with j n lying 
in the range 7„ > p is assured only for a c < a < 1, while for < a < a c , the validity of the proof of the absence 
holds for non-oscillatory normal modes with the range j„ > q. Further developments of methods to analyze non- 
oscillatory normal modes with the smaller range of j n will be necessary for claiming definitely the absence of any 
non-oscillatory unstable normal modes. However they require more complicated and sophisticated techniques to 
analyze Eq. (5.15), which also remains to be investigated in future works. 
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VI. SUMMARY AND DISCUSSION 



Now let us summarize the results obtained in this paper. We have considered a spherically symmetric system 
describing gravitational collapse of a perfect fluid with the equation of state P = ap. To treat analytically time 
evolution of the perturbations in the self -similar background (which may describe bouncing and recollapsing motion 
of a fluid), we have constructed the single wave equation governing the perturbations. As an application of the 
wave equation to the stability analysis, the eigenvalue problem for the perturbations having the time dependence 
exp {icv log(— f)} has been studied, and we have arrived at the following main conclusion: Non-oscillatory = 0) 
unstable normal modes with the growth rate y n > 1 do not exist if Eq. (5.24) holds for the self-similar backgrounds 
in the subsonic region, and some additional conditions such that H/H s < 1 become necessary for the absence of 
non-oscillatory unstable normal modes with -y„ < 1. 

In this paper we have found the absence of non-oscillatory unstable normal modes, except for ones with the small 
growth rate such as y n < q (for < a < oc c ) or -y n < p (for a c < a < 1), in the flat Friedmann background for 
< a. < 1. This supports the numerical result given by [10] for < a < 0.036. Because the attractor behavior 
of the flat Friedmann solution was also pointed out in [10], the absence of oscillatory unstable normal modes with 
nonzero f> n may be also shown in our analytical scheme. In contrast to the case of the flat Friedmann background, it 
has been numerically found that a non-oscillatory unstable normal mode can be excited in the critical gravitational 
collapse f or < a < 0.89 [13]. The consistency of the numerical result with our analytical scheme may be checked 
by showing that the condition given by Eq. (5.24) breaks down for the critical solution. 

Finally we would like to remark that the absence of non-oscillatory unstable normal modes with j n > 1 is an- 
alytically shown for the Newtonian Larson-Penston solution [17]. It is interesting to note that some mathematical 
difficulties also appear in the Newtonian analysis of the slowly growing modes with j n < 1. The result obtained by 
the normal mode analysis for the Newtonian Larson-Penston background strongly suggests that the condition (5.24) 
(for 7„ > 1) should be satisfied for its general relativistic version. 



APPENDIX A: ASYMPTOTIC BEHAVIOR OF SELF-SIMILAR SOLUTIONS NEAR THE SONIC POINT 



From Eqs. (3.10), (3.11) and (3.12) we obtain the expansion of the functions V, A' and X near the sonic point £ = £ s 

as 

OO CO 

1=1 j=0 

CO CO 

a'® = E4(£-&) i + (e-&)*E4+/(£-&V' (A2) 

;=o j=o 

CO OO 

X(0 - E X fc - {J - - ^)' T E X X ■ - > (A3) 
i=0 ;'=0 

where the variable £ is defined as £ = log(-z), and the power index % is given by 



ViX = - 



2V« 

{ (1 + a) 2 (l + 3a) Aq 2 + 2(1 + a) (1 + 5a)A' + a 2 + 6a + 1 } X x 



We also obtain 



Va(l + a)X A^ 
V X A[ 2^(l-a){3(l + a)^ + 2} 
A' x 4(l + a)A^ 

V x _ Va-(l-a) 
A' x l + (l + a)A' ' 
X x 2Xq 



(A4) 

(A5) 
(A6) 

x o = /, , ttttt: ^ 4/2 ■ r rh\ ~imn—^—T—, ■ ( A7 ) 



A' x ' l + (l + a)A' ' 

(l + a) 2 

(1 + a) 2 (3a + 1)Aq 2 + 2(1 + a) (1 + 5a)A^ + a 2 + 6a + 1 
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Moreover we can obtain the quadratic equation for A[ and find 



A ,_ A , _ (1 + 2A' ) {(1 + a)(5a - IK + 3a - 1} ± {l + (1 + a)A' } 
A 1 A 1± = ^—^ , 



where y is defined by 

y = 4(1 + a) 2 (l + 3a)A^ 2 + 4(1 + a) (1 + 3a 2 )A^, + (1 - 3a) 2 . (A9) 

Then V\ is written by 

V, = V 1±s -f (l±j^). (MO) 
Substituting Eqs. (A5), (A6), (A7) and (A8) to Eq. (A4), we arrive at the simple form 

* = -l-f- (AH) 

The solutions V, A' and X are C 1 or analytic at the sonic point when either x > 1 or A^. = (i.e., the second terms of 

Eqs. (Al), (A2) and (A3) disappear). Then, it is easy to see that V x = V/, A = A£, A x = A", X = X s and X x = X£. 

APPENDIX B: FUNCTIONS IN THE PERTURBATION EQUATIONS 

We can write the coefficient functions P, and Q, in Eqs. (4.3) and (4.4) as 

Pi = -^{l + (l+«)A'}| a V 2 A' 2 -(l + A') 2 } , (Bl) 

Pi = [l+a + X(l- a ){y 2 A ,2 +(l + A') 2 }] , (B2) 

P 3 = -2X(1 + A')A'V 2 , (B3) 

P 4 = -^{-V 2 A'(l + a + A')+cc(l+A') 2 } , (B4) 

P 5 = [l + a + x{3(l-a)V 2 A ,2 -(l + 7a) (l + A') 2 }] , (B5) 
olA' 

Ql = ITa^ 2 ' (B6) 

Q 3 = -^{V 2 A' 2 +(1 + A') 2 } , (B7) 

Q 4 = -^{-V 2 A' 2 +(l + A')(*A'-l)}, (B8) 
olA' 

Qs = TTa^ 5 • (B9) 

The coefficient function T, in Eq. (4.7) are given by 

Ti = - (P 3 Q 4 - P4Q3) f + Pi (P 3 + Qs) = -2XA' (1 + A') P X V (V ± , (BIO) 

T 2 = Pi (Q 4 + P 4 ) = 2XA' (1 + A') P x (V 2 - a) , (Bll) 

(l + a)PiP 5 A' (V±Ja) 

T3 = - (P5Q4 - P4Q5) / + Pi (ft + Qs) = A ,y±^(i + A ,) ^ ' (B12) 

where the upper and lower sign correspond to the choice of / — f+ and / = /_, respectively. In addition the function 
Y is given by 

r = (P5Q4 - P4Q5) f + (P2Q4 - P4Q2 - P1P5 - P1Q5) f-Pi (Q 4 + p 4 ) /' - ft (p 2 + Q2) . (Bi3) 
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Using Eqs. (3.10) and (3.12), we write A" and V in /' of Eq. (B13) by V and A' and X to obtain a simpler form of Y. 
Then we have 

Y= P 1 A'(1 + A')V(Y 1 X + Y 2 ) t 
^{A'V± y/u(l + A')Y 

where the functions Y\ and Yi are written as 

Yj = ±{(l + a) 2 (3a + l)A /2 + (-3a 3 + 9a 2 + 15a + 3)A'+2(-a 2 +4a + l)}A' 2 y 2 
+8aV^ (1 + A') {1 + (1 + a) A'} A'V 

T (1 + A') 2 { (1 + a) 2 (3a + 1)A' 2 + (1 + a) (a 2 + 16a + 3)A' + 2(a 2 + 6a + 1)} , 

(B15) 

Y 2 = ±(l+a) 2 {(l-a)A ,2 + (3 + a)A'+2} . (B16) 
The functions W, F and Li in Eq. (4.8) are given by 

W = -Y—^ L 2aV ' + ^ + V2K2 \ (B17) 
a \ V(V 2 -a) } }(V 2 - a) J 1 ; 

F = "v^\ — + T/' (B18) 

(V 2 - a) K 3 

U = — i i, (B19) 

a/ 

where the functions K\,K 2 , K 3 and / are written as 

JCi - (Pi + P4f)(YT[-Y'T 1 ) + (P 2 + P 4 f + P 5 f)YT 1 +YP 3 fT 3 + P 3 Y 2 , (B20) 

K 2 = (Pi + P4f)(YT^-Y'T 2 + YT 3 ) + (P 2 + P 4 f + P 5 f)YT 2 + Y 2 P 4/ (B21) 

^3 - (Pi + P4f)(yTi-rT 3 ) + (P 2 + P 4 f + P 5 f)YT 3 + Y 2 P 5r (B22) 

/ = YP 3 / T a . (B23) 

Using Eqs. (3.10), (3.11) and (3.12), we can also write the functions W, F and U by V, A' and X only. Then we have 

w = w^ + w 2 x + w 3 

v^(l+ a )VA'(l + A')(YiX + Y 2 )X' v ; 

FlX + f2 



2V^(1 + a)A' (1 + A') V {1 + (1 + a) A'} { A'V ± (1 + A') } X ' 

(B25) 



;j = Zlw± L/iX 2 + U 2 X + iJ 3 



T 4 2a(l + a)A'(l + A') V 2 {l + (l + a)A'}{A / y± v ^(l + A')}X 2 ' 

(B26) 
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where the functions W\, W2, W3, F\, F2, U\, U2, U3 and T4 are given by 

Wi = (1 + 3a) {-2(1 + a) 3 A /3 + (a -3)(l + a) 2 (3a + 2)A' 2 

+2(a - 3) (1 + «)(1 + 3a) A' + 2(a 2 - 4a - 1)} A /4 V 4 
±2Vu(l + a) (1 + A') {-(1 + a) 2 (l + 3a)A' 2 

+3(1 + a) (a 2 - 4a - 1) A 1 + 2(a 2 - 4a - 1) } A' 3 V 3 
+ (1 + a) (1 + A'f {2(1 + a) 3 (l + 3a)A' 3 - (1 + a) 2 (3a - 5)(1 + 3a)A' 2 

+4(-lla 3 - 9a 2 + 3a + 1) A' - a 3 - 33a 2 + a + l} A ,2 V 2 
±2Vk(1 + a) (1 + A') 3 { (1 + a) 2 (1 + 3a) A' 2 

+ (1 + a) (a 2 + 16a + 3) A' + 2(a 2 + 6a + 1) } A'V 

- (1 + A') 4 {2a(l + a) 3 (l + 3a)A' 3 - (1 + a) 2 (l + 3a)(l + 5a)A' 2 

-2(l + a)(l+5a)(a 2 + 6a + l)A'-(a 2 + 6a + l) 2 | , (B27) 

W 2 = -(1 + a) 2 {2(1 + a) 2 A' 3 + (6a 3 + 7a 2 + 10a + 5) A' 2 

+4(1 + a) (1 + 2a) A' - a 2 + 4a + l} A ,2 V 2 
±2V«(1 + a) (1 + A') { (1 - a) A' 2 + (3 + a) A' + 2} A'V 

- (1 + A') 2 {2a(l + a) 2 A' 3 - (a 3 + 6a 2 + 11a + 2)A' 2 

-4(3a 2 + 6a + l)A' - 2(a 2 + 6a + 1)}] , (B28) 

W 3 = (l + a) 4 {(l + a)A ,2 +2A' + l} , (B29) 
Fx = |(l + a) 2 (9a-5)A /2 + 3(l + a)(a 2 + 6a-3)A'+4(a 2 + 2a-l)| A' 3 V 3 

+Voc(l + a) (1 + A') {(1 + a) (3a - 7) A' 2 + (3a 2 + 2a - 9) A' + 2(a - 1)} A l2 V 2 

+ (1 + a) (1 + A') 2 { (1 + a) (11a + 1) A' 2 + (5a 2 + 14a + 1) A' + 4a} A'V 

+ (1 + A') 3 { (1 + a) 2 (17a + 3) A' 2 + (1 + a) (a 2 + 30a + 5) A' + 2(a 2 + 6a + 1) } , 

(B30) 

F 2 = -(1+a) 2 [{(l+3a)A' + l + a} A' 2 V + V^(l + A') 2 {(3 + a)A'+2}] , (B31) 
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Hi = ±3(1 -a) j(l + a) 2 (l + 3a)A /2 + 2(l + a)(l + 5a)A' + a 2 + 6a + l} A /4 V 5 
+6V41 - a 2 ) (1 + A') {1 + (1 + cc)A'} A' 3 V 4 
±3(a - 1) (1 + A') [ (1 + a) 2 (l + 3a)A' 3 + (1 + a)(a 2 + 12a + 3) A' 2 

+ (9a 2 + 16a + 3) A' + a 2 + 6a + l} A ,2 V 3 
+\/a(l + 7a) (1 + A') 2 { (1 + a) 2 (l + 3a) A' 3 + 8a(l + a) A' 2 

-(a 2 +3)A'-2(l + a)} A'V 2 

T2a(l + a)(l + 7a) (l + A') 3 {l + (1 + a)A'} A'V 
-V^(l + 7a) (l + A') 4 {(l + a) 2 (l + 3a)A' 2 

+2(1 + a) (1 + 5a) A' + a 2 + 6a + 1 } , (B32) 

U 2 = -(l + «) + {3(l + a)(2a 2 -a + l)A /2 +2(l + a)(3 + a)A'-3a 2 +4a + 3} A ,2 V 3 

+Voc(l + a) (1 + A') { (1 + a) A' 2 - (1 + a) A' - 2} A'V 2 

+a(l + a) + {(1 + a)A' 2 + 3(1 +a)A'+ 2} AV 

(1 + A') 2 { (1 + a) (3a 2 - 11a - 2)A' 2 

-4(1 + a)(l +6a)A'-2(4a 2 + 7a + 1)}] , (B33) 

(J 3 = T V^(l + a) 3 {V^A ,2 V± (1 + A') 2 } , (B34) 
T 4 = t2V^^XA' (1 + A') Pi . (B35) 
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